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Abstract 

in 

0^ ' We investigate the linear stability of a liydrodynamic relativistic flow of magnetized plasma in the simplest 

case where the energy density of the electromagnetic fields is much greater than the energy density of the 
matter (including the rest mass energy). This is the force- free approximation. We considered the case of light 
^ . cylindrical jet in cold and dense environment, so the jet boundary remains at rest. Continuous and discrete 

Q spectra of frequencies are investigated analytically. An infinite sequence of eigenfrequencies is found near the 
edge of Alfven continuum. Numerical calculations showed that modes having reasonable values of azimuthal 
wavenumber m and radial number n are stable and have attenuation increment 7 small. The dispersion curves 
uj = ct;(fcj|) have a minimum for fcy^ ~ 1/R {R is the jet radius ). This results in accumulation of perturbations 
inside the jet with wavelength of the order of the jet radius. The wave crests of the perturbation pattern formed 
^ . in such a way move along the jet with the velocity exceeding light speed. If one has relativistic electrons emitting 
yjf) ' synchrotron radiation inside the jet, than this pattern will be visible. This provide us with the new type of 
. superluminal source. If the jet is oriented close to the line of sight, than the observer will see knots moving 
backward to the core. 
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1 Introduction. 



Possibly the most intriguing feature of a numerous extragalactic radio sources is the existence of a narrow 
well collimated radio jets. It is these jets that are believed to be responsible for the transportation of a great 
i amount of energy from the central compact areas of the galaxies to their distant radio emitting parts. From the 
H I observations of the superluminal motions of bright knots along the jets one must conclude for the velocity of 
the flow to be relativistic with the Lorentz factor 7 being of the order of 5 to 10. One of the important problem 
in the physics of extragalactic jets is their stability over the large distances. Up to now there have been many 
works in which the stability of the jets is investigated under the different assumptions for the velocity of the 
matter in a jet and the influence of the magnetic fleld on the flow dynamics. Turland & Scheuer (1976) and 
Blandford & Pringle (1976) were the first who considered the Kelvin-Helmholtz instability of a plane relativistic 
hydrodynamic flow in vortex sheet approximation and applied the results to extragalactic jets. The perturbed 
modes are considered both for cylindrical geometry of a jet and for the plane boundary between the jet and the 
outside media if the wavelength is small enough (Hardee 1979, 1987; Birkinshaw 1984; Payne & Cohn 1985). 

When it had been recognized that the magnetic fleld may be dynamically important for jet confinement, 
the problem of magnetohydrodynamics (MHD) stability of cylindrical flows naturally arose. This facet of the 
problem goes back to the investigations of the stability of magnetostatic plasma equilibria. Appert, Gruber 
& Vaclavik (1974) have pointed to the existence of Alfven and slow wave continua in the case of the sheared 
magnetic field and (or) nommiform plasma density. Cohn (1983) found the uniform jet with the "top hat" 
velocity profile confined by the poloidal current flowing on the boundary, unstable for all modes. He has also 
investigated relative importance for instability fundamental and reflection modes (these modes differ by a radial 
wavenumber n, fundamental mode is one having n — 0, reflection modes have n > 0). The general system of 
the two flrst order differential equations on the radial Lagrangian displacement S,r of the fluid element and the 
perturbation of the total pressure Pi , which describes the propagation of small disturbances along the cylindrical 
nonrelativistic MHD flow, were derived by Bondeson, lacono & Bhattacharjee (1987) (thereafter referred to as 
BIB). Interesting information concerning the existence and location of point eigenvalues in the complex lu plane 
was obtained in BIB from a boundary layer analysis near the Suydam surface kB — 0. Recently, Torricelli- 
Ciamponi & Petrini (1990) and Appl & Camenzind (1992) investigated numerically instabilities of nonrelativistic 
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MHD cylindrical flows with nonzero poloidal magnetic field and poloidal current distributed over the jet cross 
section. 

Next in order of complexity to relativistic hydrodynamic and nonrelativistic MHD approaches is relativistic 

MHD case, which is the generalization and unification both of them. When considering relativistic models with 
a significant magnetic field it is necessary to involve the electric force and, for non-stationary case, to add the 
displacement current. Indeed, at the speed of the hydrodynamic flow u ~ c the electric field in plasma with high 
conductivity is of the order of the magnetic one E = — v x B/c. The charge density g is equal to V-E/47r ~ E/L 
and the current density j in a stationary flow or for small nonstationarity is j = c • (V x B)/47r ~ cB/L. We 
see that the ratio of the electric force, £>E , acting on the unit volume, to the magnetic force, j x B/c, is of 
the order of E'^ / ~ f^/c^ ~ 1 for the relativistic case. So it would be, probably, better to use the word 
electro-magnetohydrodynamics instead of magnetohydrodynamics for relativistic flows. Many extragalactic 
and, possibly, stellar jets (see recent discovery of superluminal motions by Mirabel & Rodriguez, 1994) appear 
to be relativistic and magnetized, so there are numerous papers devoted to the construction of selfconsistent 
description of such flows (Camenzind 1987; Lovelace, Wang & Sulkanen 1987; Li, Chiueh & Begelman 1992; 
Contopoulos 1994). It is of interest to examine the stability of the solutions obtained. 

We started this investigation with the simplest case where the energy density of the electromagnetic fields 
is much greater than the energy density of the matter (including the rest mass energy). This is the force- free 
approximation. This case is reversal to the pure hydrodynamics one when the electromagnetic field is absent. 
Using force-free approximation one can hope to take into account the influence of the electrodynamic effects 
on the stability and simultaneously obtain substantially simplified problem. The terms in the momentum 
equation which are proportional to the mass and the pressure of liquid, are therefore small compared to the 
electromagnetic force pE + j x B/c, so it is possible to set pE + j x B/c = (we use imits where c = 1). 
The force-free approximation together with the ideal hydrodynamics approximation (which means an infinite 
conductivity of plasma and consequently the absence of the electric field in the frame moving with the clement 
of the medium) can be applied to the neighbourhood of the massive black hole, which is thought of as the 
central engine of active galactic nuclei. Such an approach was developed by Blandford & Znajek (1977) and 
Macdonald (1984) [ see also chapter IV of the book "Black Hole: The Membrane Paradigm" by Thome, Price 
& Macdonald (1986) and chapter VII of the book "Physics of Black Hole" by Novikov k Frolov (1986) and 
references therein]. 

In our previous paper (Istomin & Pariev 1994, thereafter referred to as IP) we consider the stability of a 
force-free jet with respect to axisymmetric perturbation. Present paper is the continuation of IP and deals with 
nonaxisymmetric perturbation. As in IP the equilibrium configuration of the jet has cylindrical symmetry. This 
means that all quantities describing the jet depend on the distance from the jet axis r and do not depend on 
the coordinate along the jet z and rotational angle (j). The boundary of the jet has the shape of a cylinder. We 
suggest that the jet propagates in the medium which density is greater than that of the jet but temperature 
and pressure are small, so the condition of impenetrability is fulfilled and the boundary is at rest. The poloidal 
magnetic field B^ is assumed to be uniform and parallel to the jet axis. The fluid moves along spirals because 
of the radial electric fleld. It has been proved analytically (IP) that under such conditions the relativistic 
flow is stable for axisymmetric modes (to = 0, where to is azimuthal wavenumber). In section 2 we derive 
equations governing the problem, describe the procedure of finding eigenfrequencies based upon the Laplace 
transformation method, and discuss asymptotic behaviour of perturbations in a long time after initial excitation. 
In section 3 we present the results of numerical calculations. In section 4 we perform boundary layer analysis 
of our equation near the point va in the complex plane r, which is the point of the coincidence of two Alfven 
resonant points. Possible astrophysical implications of the results obtained is discussed in section 5. 



2 Stability problem 

Let us consider a fiow of liquid in a force-free cylindrical jet. We use the units in which c = 1. The condition 
of an ideal flow is 

E = -vxB, (1) 
where v is the plasma velocity. The force-free approximation is guided by the relation 

£.E+jxB = 0. (2) 

First we will review the stationary configuration of the jet described in IP. In this case V x E = and the 
velocity v can be written as 

v = KB + n^re^. (3) 
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Here and below r, z and </) are cylindrical coordinates, e^., and are the unit vectors in the cylindrical 
coordinate frame, K = K{r), ^l^ = U^{r). Then 

E = -n^r{e^ X B) (4) 

and can be treated as the angular rotation velocity of magnetic field lines (Thorne et al. 1986). In the 
cylindrical configuration B = Bz(r)ez + B^{r)eff,, so using Maxwell equations, we obtain from (|^) 

E(V • E) - B X (V X B) = 0. (5) 

According to (^) the only non-zero component of equation is r-component. This implies 

n^BM^^r'Bz) = + lB^±irB^). (6) 

ar ar r ar 

This equation governs the force balance in radial direction and defines all possible solutions for the force-free 
electromagnetic fields in cylindrical configuration of the magnetic tubes. General solution of this equation was 
described in IP. For uniform poloidal magnetic field B^ = const, we have 

B^^±Q^rB,. (7) 

If B^ is constant the stationary magnetic field structure is entirely determined by the function fl^{r). 

The case where the total charge of the jet is equal to zero is probably the most natural. If the jet has 
a charge not equal to zero, the electric field penetrates into the surrounding medium. This results in charge 
motion in the plasma and in a decrease of the charge of the jet. To avoid the problem of closing current loop 
somewhere outside the jet it is naturally to demand the total poloidal current through the jet to be equal to 
zero. These two requirements lead to fl^{R) = 0, where r = i? is the jet boundary (IP). The equilibrium 
stationary configuration of the jet is shown in Fig. 1, which is extracted from IP. We reproduce it here for 
illustration purpose. 

2.1 Basic equations 

We perform linear stability analysis using common method of small perturbations. In the subsequent formulae, 
the values referring to nonperturbed solution will be denoted by the subscript '0', while ones referring to 
perturbation - by the subscript '1'. We will consider only nonaxysimmetric perturbations, so throughout the 
rest of the paper it is assumed m ^ unless specified directly. 

After removing the quantities g, E, j from the initial system of Maxwell equations and equations (|l|)-(^) 
there remains only three resultant equations involving B and v only: 

V • B = 0, (8) 

— - V X (v X B), (9) 
(v X B)V • (v X B) - B X (V X B) - B X 

^(v X B) - 0. (10) 

The quantities B and v can be represented as B = Bq + Bi , v = Vq + vi. We consider the perturbations in 
the form 

Bi = bi(r) • exp(— iwt + ikz + im(f)), 

Vi = ai(r) • exp(— iwt + ikz + im(f)), (11) 



where m is an integer. Substitution this expressions into linearized set of equations (pD-QlOD and removal the 
components of the perturbed quantities with subscript '1' in favour of Bri lead us to the following second-order 
ordinary differential equation of the variable Bri, where the prime denotes differentiation with respect to r: 



„ „ , I'l 2m2 1 dn'' 1 , 

Bri + Bri 5 2 71 m ~ "T"; OF + 



2 ,2 

uj - k — 



p d f 1 



^ ' dr \{k + mn^)^ dr 

{AiBr'i+A2Bri)=Q. 
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2mnP -nP^r^{uj + k) 



(12) 
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The quantities Ai andA2 are 

A, = -2nP\{u + k)+ '^-^i±JL-{m - mn^\^ - 2kr^n^) , (13) 

dr k + m\t^ 



Ay = 2m 



dn^\^ n^r^{uj + k)-m dfl^ 1 (2m 1 



dr J k + mQ^ dr k + mQ^ \ r uj^ ^ k'^ — m^ / r' 



— {3mn^ -LO + k)-n^ r^LO^ - P) ln''iuj + k) + — 



+ -{cj-k)+ (14) 



I X ?T? 1 

2n^ rk{uj + k) + rm{LO + k)) + —VL^ + —{uj - k) + 2,^1^ {lo + k) 



2 1 ^2 1 / ^2 



771 1 

r^LO + k r^ oj^ - k^ - m^ /r^ \r^{uj + k) ^ ' 

This equation has been obtained in IP. When deriving this equation it was accepted that B^q = const and 
the relation (|^) with the sign '+' was used. The choice of the sign in formula (|^) depends on the choice of 
the direction of the poloidal component of the flow in the jet. This keeps equation ( p^ with definitions ( p^ ) 
and ( p^ unchanged under reversal of the signs of B^^ and k. For definiteness we adopt the sign "+" in (^ and 
consider arbitrary values of k. All other perturbed physical quantities can be readily expressed through Bri- 
The expressions for some of them are given in the Appendix A. Briir) must fulfill the boundary conditions for 
r = and r ~ R. For r = 0, B^i must be regular. If |m| ^ 1 this condition can be strengthened to become 
Bri\r=o = 0. The boundary condition for r ~ R must be derived from the rigidity of the jet wall Vri\r=F!. — 
which has been assumed in the present investigation. It gives Bri\r=R = 0. Thus, in order to find perturbed 
modes, we have to solve the edge problem for equation ( p^ ) with the boundary conditions 

Br ilr^o is regular, 

As a result of calculation parameter K{r), which determines the component of Vq parallel to Bq, drops out 
from the equation ([l^). Therefore, the results of our investigation of stability do not depend on the values 
and profiles of the longitudinal velocity. The physical reason for that is the neglecting in the force-balance 
equation (^), which is essentially the Eueler momentum equation, the terms describing the contribution from 
the inertia of the mass fiow and the pressure. The results of the stability investigation do not depend also on 
the conductivity of the outside medium, since there are no perturbations in that region. 

The primer goal of our investigation is to answer the question whether the jet is stable or not. So it is more 
convenient to use temporal approach for investigating stability, i.e. seek for complex values of u; for real k. 
The existence of only one eigenvalue lu with positive imaginary part would mean that the jet configuration is 
generally unstable. This is not the case for spatial approach, when the search is carried out for complex values 
of k for real lo. In this case the sign of imaginary part of k is irrelevant for stability and more complicated 
analysis should be done to answer the question of stability (Lifshitz & Pitaevskii 1979). We adhere here to 
temporal approach. 

The radial eigenvalue problem ( p^ ) has regular singularities in the following 4 cases: 

1. k + mil^ = or kBo = 0. This is the resonant surface r — re- For nonrelativistic MHD consideration 
there exist local Suydam modes in the vicinity of Tc (BIB). However this is not the case for the force- free 
approximation (we see it below). Expanding equation (^2|) to lowest order in x = r — Tc, substituting 
Bri (X x'^ we find the indicial equation {v — l)(i^ — 2) = 0. This implies that the solution near r = rc 
is the linear combination of two linear independent solutions, Wi{x) and W2{x), with some constants 
ki and ^2 B^i = kiWi{x) + k2W2{x). The first terms in expansion of these solutions are Wi{x) = 
aix^ + a2X^ + a^x^ + . . ., 1^2(2;) — AWi{x) \ogx + bix + 622;^ + b^x^ + . . . It can be shown directly 
by inserting the expression for W2{x) into equation (p^ ) that A ^ 0. Therefore, this singular point does 
not lead to nonanalytic solutions in it's neighbourhood and is only apparent. 

2. ti)^ — fc^ — m^/r^ = 0. This point can be interpreted as the resonance of the perturbation, having k parallel 
to magnetic flux surfaces of Bq, with the fast magnetosonic wave which in the force- free approximation 
always has the speed equal to the speed of light. Indicial equation is — 2) = 0. Similar to the type 1 
singular point it can be shown that the coefficient of log a; in the expansion of Bri near resonance vanishes, 
therefore this singular point is only apparent too. 
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3. k-Lu + 2mn^ ~ r^{uj + fc) = 0. This is the resonance of the perturbation with the relativistic Alfven 
wave r = TA- When l.h.s. has a simple zero, the indicial equation will he = Q and, in general, -B^i 
has logarithmic singularity there. The expansion of other physical quantities near r = is given in the 
Appendix A. 

4. r = 0. This is the singularity produced by the coordinate origin. Two linear independent solutions near 
r — Q are Bri (x rl"*'"^ and Bri oc r~l™l~^, where m ^ 0. In order to meet the boundary condition ( p^ ) 
one should choose the former solution only. 

One can solve edge problem for equation ( p^ ) directly as it was done in IP for the case m — 0. However, 
the case m ^ is more complex, particularly, because of the existence of the singularities in the coefficient of 
equation mentioned above, two of them being only apparent. The numerical treatment of equation 
can be greatly simplified and integration through 1-st and 2-nd type singular points can be avoided if wc rewrite 
equation (|l^) as a system of the two first order differential equations in terms of the radial displacement and 
the disturbance of the total pressure Pi instead of the radial component of the magnetic field perturbation Bri- 
We need transformation into variables and Pi to consider initial value problem in subsequent subsection and 
also in order to possible further generalization of the problem to include inertia of plasma. For the Lagrangian 
displacement (i.e. the displacement of a fluid element moving with the equilibrium flow) ^(r, t) one has well 
known expression relating it to the Eulerian disturbance of the velocity field vi 



8$ 

at 



vi + (|V)vo - (voV)|. 



(16) 



The r-component of this relation is (bearing in mind the representation ( |ll| 

m 

- iuS,r = Vri - ikVzQ^r - 'i — 4> o^r ■ 

The r-component of the linearized induction equation (^ is 

I m \ ( m \ 

\^-UJ + kv^o + — W0oj ^rl = [kB^Q + —B^i^j Vri- 



(17) 



(18) 



Combining equations (|lj) and ( p8| ) one can readily obtain that Bri = i {kB^a + ra/rB^^ which in the 
force-free case with Bzq — const transforms to 



Bri = iS,rBzoik + mQ^). 



(19) 



We call the value Pq = (Bq — E§) as a total pressure because this is the diagonal part of the electromagnetic 
stress tensor tTy — gt- (ij2 _ jtj2^ . _ _L (^Bi^Bj^ — Eii^Ej^. Next the perturbation of the total pressure will 
be ^ 

Pi — ^ {BzqBzI + B^^Bcf,-^ — EroEri) , (20) 



where the radial component of the perturbation of the electric field is Eri 



-{BzqV^^ — B^^^Vzl + V^qBz 



VZ0B4, 



) and radially directed equilibrium electric field is ErQ — VzqB^^ — v^^Bzq — —^ rBzQ- Using formulas 



(A1)-(A3) from the Appendix A, which express B^f,^, Bzi and Eri by means of Bri and dBri/dr, we find for 
Pi expression by means of Bri and dBri/dr 



iirPi = iB, 



S 



AdBri „ f dn^ A If , 2/ 

■ Bri<m——— + — (uj - k~n^ r\u} + k) 
dr r ^ rl' \ 



F dr 



(21) 



where 



F = k + mn' 



A^k-, 



1 2 111 

k — m /r , 

^ + 2mn^ - n^^r^iuj + k). 



(22) 
(23) 

(24) 



Those values of r, which make F, S or A equal to zero, are 1-st, 2-nd or 3-d type singular points respectively. 
From ( |2l] ) using ( p^ ) dBri/dr is expressed by means of Pi and ^r 



dBri 

dr 



AttPi 



iFS 



(uj + k)BzoA 



^oj — fc — Q^^r^{LU + k)j + mr 



dn^ 

dr 



BzoCr- 



(25) 
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Differentiation of the equation ( |19[ ) and substitution the result for dBri/dr in ( P5[ ) lead after some cancellation 
to the first equation from the couple desired 



dr 



S 



B,l A{uj + k) 



k) 



(26) 



To obtain the second equation, which relates dPi/dr to the linear combination of Pi and ^i, we differentiate 
equation on r, substitute for appeared in r.h.s. of ( pl| ) d^Bri/dr^ it's value from the second order differ- 
ential equation ([l^), and remove dBri/dr and Bri in favour of Pi and by means of (|25|), (|l9|). Numerous 
cancellations occurred when this procedure was carrying on. Finally, we get the following equation 



47r 



dPi 2n^ 



dr 



A 



(n^r{uj + 47rPi - (w + fc) ( ^ 



A 



(27) 



We introduce for convenience the dimensionless pressure disturbance p.^, 
tions (|2^) and ( p7| ) can be rewritten as follows 

r dr 

A^^Csir-Cip,, 
dr 



AnPi/B^l. Than pair of equa- 



where 



Ci 

C3 



m'^/r^ 



uj + k ' 
= -(cj + fc)(A2 -4f7^^). 



(28) 

(29) 
(30) 
(31) 



The system (|2^) has the same form as derived by Appert, Gruber & Vaclavik (1974) for nonrelativistic 
MHD stability investigation of the static plasma cylinder. The only difference is in the coefficients A, Ci, C2 
and C3. The presence of the equilibrium flow does not influence the form of the set of the equations, but only 
modify the coefficients (e.g. BIB). Now we see that for relativistic force-free flows this result remains to be 
true. It is easier to integrate numerically set of equations (^ ) than the second order equation ( [l^ ) and in our 
numerical computations we actually integrated (|2^). First, one should note that no derivatives of fl^ (r) enter 
into the (|28|). Second, which is more important, the only singularity in ( p8| ) is A = 0. The values of r which 
makes S" = or F = are regular points of the system (28). Because of the relation ([l^ ) it is clear that they 
will be regular points of the second order equation ( [l^ ) on Bri as well. This result is just that we have obtained 
above by considering explicit expansion Bri the neighbourhood of singular points. 

Equations (28) can be readily converted into one second order, so called Hain-Liist (Hain & Liist, 1958) 
equation on r^r 

'd_ flCi\ C5-C2CA 
dr\rC2j rC2 



dr \C2 r dr 



r^r 



= 0, 



(32) 



where we made notice of a certain factorization 



Ci — C2C3 = A{C5 — C2C4). 



The C4 and C5 are 



-A{uj + k), 



C5 = -4f7^ (tu + k). 



(33) 

(34) 
(35) 



Equations (p2[ ) and (33) are identical in form to that derived in BIB for nonrelativistic MHD. As well as in 
equation (|12[) in Bri the zero C2 = in ( [32[) is only apparent singular point, S,r is regular at this point. If the 
factorization (|3^) did not occur, equation ([32[) would have essential singularities when A has a zero of quadratic 
or higher order. 
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2.2 Laplace transformation method 



Now edge problem (jlj) , ( |15[ ) is reformulated for system ( |28| ) or equation (|3^) with evident boundary conditions: 
is finite at r = and = at r = i?. Thereafter we shall use dimensionless values r' = r/R, lu' = cuR, 
n'' = n'^R, k' = kR and the prime will be omitted. Thus, r = 1 will correspond to the jet boundary. When 
integrating equations (|2|) the problem arises how to treat the singularity ^ = 0. To answer this question it 
is necessary to remember that we use temporal approach, i.e. solving initial value problem. We seek for the 
solution ^r(i,r) for all t > having given the initial conditions, say, vi(0, r) and Bi(0,r). Because of the all 
equations governing the problem is linear on perturbations and the explicit dependence of this equations from 
4> and z is absent, it is enough to consider the behaviour of only one Fourier component of all disturbances with 
respect to (j) and z. As well as in previous text we shall take the perturbation and initial conditions having 
the form (,rkm oc exp{ikz + iuKp) and shall omit the subscripts m and k belonging to all the disturbances and 
initial conditions. Instead of doing Fourier transformation with respect to t it is useful to apply the technique 
of Laplace transformation. We introduce function ^ruii'r) determined by 



^ruir) = J Ut,r)e dt. 



Than the reverse Laplace transformation is 

Ut,r)= J e™Me-^-*^, (36) 

where the integration is performed along the line in the upper half of the complex plane uj, which is parallel to the 
real axis and goes above all irregular points of ■ For the sake of resemblance to the Fourier transformation we 
use here the frequency w related to the usual parameter p of the Laplace transformation by w = ip. The Laplace 
transformation of other perturbed values is defined in the same manner. Carrying out the Laplace transformation 
of the equation (^6|) and the source equations of the problem eliminating all disturbances in favour 

of £,ruj, vi(0,r) and Bi(0,r) we get finally the equation ( ^2|) with some r.h.s., which depends on the initial 
conditions. 

Note, that we need to add to vi(0,r) and Bi(0,r) the initial condition for the displacement ^(0,r), when 
processing with the (|lq). Contrary to that, if we did not introduce displacement and did restrict ourself to 



Laplace analogue of equation (12) in Bri, the given vi(0, r) and Bi(0, r) would perfectly determine the unique 



solution for vi(i) and Bi(t). Thus, instead of (|28|), the equations for and p^,^^ will be 

r dr 

- Cs^r^ - C,p,^ + D2, (37) 

dr 



where Di and D2 are lineary dependent on the initial conditions and do not contain any denominators except 
u! + k and Doppler shifted frequency —a — uj — kvQ^ vq^ — uj — mfl^ — KB^^ik + mO,^). Further, we derive 

second order differential equation on rS^ruj similar to ( |3^ ) but with some nonzero r.h.s. Reduced to the normal 
form it will be 



^ , rC2 d fl A\ d . ^ . ^ 
d^^^''^-^^—Tr[-rC-2)Tr^'^-^-''^^ 



rC2 d (}_Ci\ C5 - C2C4 
A dr\rC2)^ A 



-C.d (D,\^r^^^^^^^^^^ (38) 



A dr\C2 A^ 



To solve edge problem for (^8|) with boundary conditions r^ruj\r=Q oc r"'™' , r^rw|r=i = we do the following. 
First, integrate equation ( ^ ) (or, which is equivalently, system (^^) from r = to r = 1 with the initial 
condition r^rcj|r=o = r"'™'. Thus, we obtain the solution gQ{r,iL!,k,m) of nonuniform equation (pq), which, in 



general, does not satisfy the boundary condition at r = 1. Then we integrate uniform equation (32) starting 
from r = with the same initial condition r^rcj|r=o = r^'™'. As a result we obtain some solution gi{r,uj,k,m) 
of (^). If gi is equal to at r = 1, than this will be an eigenfunction of (^) with the frequency lu being 
the eigenfrequency of the problem ujnm{k). The integer number n counts different eigenfrequences with the 
same m and k. We assume the muchness of eigenfrequences to be countable. This corresponds to the discret 
spectrum. Usually there exists an infinite sequence of ujnm for given m and k. This sequence has either finite 
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or infinite accumulation point in w-plane or both (see numerical results below and section 3). The numbers n 
can be chosen so that the higher n the more oscillations has the corresponding radial eigenfunction. The least 
oscillating mode is referred to as a fundamental one in literature and the number n = is usually ascribed to this 
mode (Birkinshaw 1984, Payne & Cohn 1985, Appl & Canienzind 1992). All other modes are known as reflection. 
In the case u! = LUnm{k) the edge problem for the equation (|3^ ) has no solutions unless go{l,uj,k,m) = 0. If 
go(l,uJnm, k,m) — than the number of the solutions would be infinite, because of all functions go + cgi with 
arbitrary constant c would be equal to at r = 1. Then for lo not being the eigenfrequency the edge problem 
for ( |38| ) has unique solution 

^™ = go{r,uJ,k,m) ; -gi(r,uj,k,m). (39 

gi{l,UJ,k,m) 



To find out ^r(^, ?") one needs to make inverse Laplace transformation (36). The asymptotic behaviour of 
^r(i, r) for t +00 is not known a priory, therefore one can be sure that the definition of S^ruj is meaningful only 
for uj having Imuj +oo. It means that the integration in ( |3^ ) must be done along the contour having a +oo 
(long dashed line in Fig. 2). In this region of complex w-plane S,ru! must be an analytical function in lu and 
it's values can be found from ( ^9| ) by applying the procedure described above. Then, the contour integral ( |36| ) 
provides us with the full solution of the initial value problem. To answer the question of stability it is enough 
to find only the asymptotic behaviour of ^rit,r) for t — s- +oo. For this purpose it is necessary to continue 
analytically the expression ( |39| ) into the whole complex w-plane. Let us fixed some value of r = (0 < < 1) 
and consider the singularities of as a function of complex uj. We assume that the initial conditions vi(r), 
Bi(r), the angular rotation velocity of magnetic field lines fl^ , and, therefore, Di{r,uj), D2{r,uj) are an entire 
functions in the whole complex r- plane. In w-plane the only singularities of Di and D2 are the poles cj = — fc 
and to — mVl^ + KBoz{k + mfl^), while the coefficients in l.h.s. of the equation ( p8| ) remain to be regular at 
these points. 

First of all, will be irregular if go{r^,,Lu) or gi{r^,,Lu) have singularities. It is seen from ( |37| ) that this can 
happens when A{r^,Lu) — 0, i.e. when Alfven resonant point coincides with chosen r*, when lo = — fc, and when 
Doppler shifted frequency a is equal to at the point r*. Notice, that from look at the equation ( |38| ) one might 
believe C2(r, w) = to be the singularityof go- However, it does not because of this is not the singularity of the 
original system ( |37| ) . From expression (^) for A one can readily conclude that for any complex uj the solution 
ryi(w) of A{rA,uj) — must be complex too. Therefore, ta can be real only for real values of lu, and gi, go are 
regular in the whole upper complex uj half-plane for any real r. By expanding (^sj) near the point r — rA{Ld) 
one can deduce the following behaviour of go and gi in the neighbourhood of the Alfven resonance 

go = cio log X + C20+ C30 log^ X + C40X + C5ox^ + ... , 

gi = cii log a; + C21 + C31X + C4,ix^ + ... , (40) 

where x = {r — rA)/rA, and Cij(w) are constants with respect to r. We see that r = is the logarithmic type 
branch point of gi, go , and consequently, ^ruj, considered them as the functions in the complex r-plane. On the 
other hand, for fixed r — r^, and uj close to the Alfven resonant frequency a;^(r*) (that is A(r,, Wyi(r*)) = 0) 
the value x can be expressed as 

^=!i(^_^^) + !i ( J4_!:k) {u:-UAf + ..., (41) 

where 



dvA 
cLlo 



(Pta 



dLJ^ 



Inserting (^) into expressions ( [40[ ) one can readily see that uj — LUAir*) is the logarithmic type branch point 
for go, gi, and, therefore, for ^rij(?'*,w). From expression ( p^ ) for A we find 

u..(n)^^+^"-""^^-)-""'^--)--^ (42) 

1 + I7f^'(n)r2 

This equation indicates that for any real function S^rio has unique branch point due to the Alfven resonance 
U!A {r* ) and this point lies on the real axis in the complex w-plane. At uj = —k function £^ruj lias a pole. The reason 
for this pole is the resonance of the perturbation with the Alfven wave, propagating in the negative z-direction. 
In the case of uniform poloidal magnetic field B^o = const this wave has always the velocity equal to the light 
velocity irrelevant to the value of and to the curling angle of the magnetic field lines Vl^r^, (see Appendix B). 
In the case of B^o 7^ const Alfven waves propagating in both directions of z-axis become equal, both their 
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velocities do depend on and the second Alfven resonance point aj^(r*) appears. It's analytical properties are 
the same as for the first Alfven resonance point: it is a logarithmic type branch point for ffii and, therefore, 
for ^rw(''*,w), with the same expansion of go, gi near it as the expansion (^0|). It can be shown that in the 
general case of nonuniform poloidal magnetic field Bzq both Alfven resonant frequencies are always real 

for real values of r* (see Appendix B). The singularity a{rf,uj) = 0, or a; = w/(r*) = mf]^ + KBoz{k + mft^) 
is also real for any real r■^, . It corresponds to the convection of the displacement ^ with the fluid velocity vq . 
oj = uj f{r^,) is a. logarithmic type branch point for go and £,rujir*,uj)- 

The second reason for irregularity of £,rui in the complex w-plane arises when = of r^(ti;) = 1, 

i.e. when Alfven resonance surface is located on the jet axis or on the jet boundary. In the first case it is 
impossible to fulfil the condition r^^uj |r=oc< r"'™' for any solution of uniform equation (|3^). When r^(cj) = 
and ^^\r=o 7^ the behaviour of the two linear independent solutions near r — are 



Therefore, in this case any solution of (^2|) will infinitely oscillate and grow when approaching r — 0. In the 
second case we have from ( ^o|) that in the vicinity of r = 1 go/gi oc \ogx, where x is the same as in (|40|). So 
solution (|9|) can not being constructed. The points 

ujAiO) = k + 2mn''{0) (43) 

and 

^ fc + 2n.l^^(l) -m^^l) 

are branch type singular points of i^rw 

At last, ^rLj has also singularity whenever lu coincides with the eigenfrequency LOnm{k). If uJnm{k) is not equal 
to liJa{^) then 171(1, w, fc, m) can be expanded in lo into Taylor series in the vicinity of the point lo = ujnm{k) 



3i(l,w,fc,m) = — 



(w - LUnra) + 



{UJ - Unraf + • ■ ■ ■ (45) 



Now we see from general solution ( p9[ ) that the function has the pole at u; = u)nm{k) unless 5o(li "t-) = 
0. If .go(l,ix') has at u; = LOnm the zero in lo of the order not less than that for 31(1,0;), then r^rw is regular at 
LO — LOnm- The casc of absence of the pole at eigenfrequency point takes place only for specially chosen initial 
conditions such that the eigen mode corresponding to LOnm is not excited. For arbitrary Vi(0,r) and Bi(0,r) 
function £_r<^ has poles at every eigenfrequency LOnm{k)- 



2.3 Analytical properties of 

Let us turn to the analytical continuation of S,ru){f*) from the region Imtj — > +00 down to lower values of Im. lo. 
We continue to consider the Laplace transformed displacement ^^i.^ at some fixed surface r — inside the jet. 
From the consideration above it follows that in the upper complex lo half-plane {i"* ) may have only poles LOnm 



and singular points ryi(w), ff{uj) of the system (37) never become real valued. Hence, the analytical continuation 



of Ct'uj(''*) into the whole upper lo half-plane (except poles lo = LOnm) is provided by the formula (t59|) without 



any changes, i.e. when obtaining gi and go integration of the system (37) should be performed along the real 
axis in the complex r-plane from r = to r = 1. This is, however, not the case for further continuation of 
C™(f'*) into the lower lo half-plane. For Ithlo — the ryi(aj) cuts the real r axis and moves for Imw < into the 
opposite complex r half-plane. Note, that form expression ( p2| ) it follows that the real valued frequency always 
exists, the crossing point lying between and 1. We want to stress that in generally not all the points rA{io) 
cut the real r axis when Ivcilo becomes 0, some of them may always lie in the upper r half-plane, another — 
in the lower half-plane. But if it occurs that ryi(tj) is real than the frequency lo, at which it does, is real too. 
To perform the analytic continuation we must deform therefore the path of integration in r of the system (|37| ) 
into a contour extending into the complex r-plane and going around the point r^(a;). If r^(w) cuts the real r 
axis beyond the interval (0, 1) then integration path is not deformed. We also need to go around the points 
r/(w) when calculating g^. After determining the integration path in such a way we calculate ^™(^*) according 
to formul a (^9| ), where now the functions go and gi are the integrals of the (^^ (or, equivalently, second order 
equation ([38[)) along deformed path. 

The procedure of analytical continuation is illustrated by Fig. 2. Initially integration in the formula of 
reverse Laplace transformation ( p6| ) is performed along long dashed line. Function ^^i^ is continued along the 
paths indicated by short dashed lines. Letters A, B, C and D denote the paths corresponding to 4 possible cases 
of the location of their intersection points with real r axis. Corresponding locations of Alfven resonant points 
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rA{oj) and their trajectories when changing cu are shown on Fig. 3 by the paths A, B, C and D respectively. In 
the cases A and D cuts the real r axis of the complex r-plane beyond the interval (0, 1) and the contour 

of integration of equations ( |3^ ) remains to be the interval of the real axis (0, 1). In the case B rA{i^) cuts 
the real r axis between and and drags the contour of integration. In the case C cutting point lies on the 
interval (r», 1) and the contour of integration is dragged from the other side of the point r*, at which we have 
to determine 

The procedure of analytical continuation described does not give unique results for the lower uj half-plane. 
The value of £,ru/{r^) depends on the path in w-plane along which the continuation is performed. Indeed, we 



may construct the solution given by the formula (39) not only for real values r but also for complex one 



For this purpose when obtaining go{r) and 51(7') it is necessary to perform integration along the contour in the 
complex r-plane, which connects three points: 0, the point r^, at which we are interested in the ^ru, and 1. By 
this way we obtain an analytical continuation of ^rLjif) into the complex r-plane. The only singularities of such 
defined ^ruif) are logarithmic type branch points r = rA{Lo) and r = rf{uj). Hence, for analytical continuation 
of to be unique, one should do branch cuts attached to the points rA{(jj) and rf(u)). Once chosen when 

luiLo — !■ 0, these branch cuts are drawn with the points ryi(a;), rf{uj) when analytical continuation of ^^(r) is 
performed in lu. Contour of integration of equations ( p7| ) in r-plane must never cross them. These branch cuts 
are shown on the Fig. 3. If we find only gi (which is enough for determining eigenfrequencies ujnm), than we 
shall not pay attention to the points r/(ijj), because they are only in the r.h.s. of equation (p^). It is seen that 
the point r^, falls on different sides of the branch cut when cu is changed along the paths B and C. Hence, if 
we continue ^rw(''*) along the path C in w-plane, we will obtain one value, while along the path B the other. 
Continuations along paths B and A give different results also. The reason for this is that in the case B the 
contour of integration goes around the singular point r — r^(w), while in the case A it does not, therefore, we 
obtain different values of 50 and gi at the point r = 1 needed to be known in formula (^9|). Because of the same 
reason the results of the continuation of ^^c^ along the paths D and C are also different (but they are the same 
for paths D and A). Thus, for the sake of ^ruj being well determined in the lower lo half-plane we need to cut it 
and keep the continuation paths from crossing the cuts. The cutting and choosing the continuation paths may 
be done in different ways. In our computations we choose the paths to be straight vertical lines in the w-plane, 
i.e. along the path Imuj runs from +00 to desired value, while Hecu remains fixed. The branch cuts are also 
straight vertical lines going down to infinite negative Imw. Our choice is refiected in Fig. 2. 

Such determined may have another singularities in the lower lo half-plane. First of all, there can exist 
the poles at the discrete eigenfrequences of the stability problem to — uJnm with linLOnm ^ 0. Secondly, the 
singularity may arise at some lo = coc when the two points r^(w) merge together. In this case contour of 
integration in r-plane becomes clutched between them and the integration procedure will be undetermined. If 
we consider two paths in w-plane with slightly different Rew such that Rewi < Rewc < Rew2 (paths C and E 
on Fig. 2) and plot corresponding trajectories of rA{i-o) in r-plane (see Figs. 4 (c) and (e) for illustration) we 
find that certain "reconnection" of trajectories occurs near the point r^(ajc)- For Imtu < Imojc the contour of 
integration catch on different r^(w) when they are moving away from merging point rA{LOc) each on it's own 
trajectory. Therefore, the values of 30(1)1 .91(1)1 C™('"*) for near loi and UJ2 will in general strongly differ from 
each other, and cUc will be a singular point of ^rui with the corresponding cut to be attached to (see Fig. 2). We 
shall see below that ujc is an accumulation point of the poles. It is necessary to stress that the singularity of 
^ruj arises only if the contour of integration goes between the merging points r^(ti;). If it does not go around no 
one of them or bypass both than lo — lOc is a regular point. In particular, all merging points having Imcoc > 
are regular. The singularity may arise only for those having Im ujc < 0. 

2.4 Discrete spectrum and Alfven continuum 

Let us now consider the asymptotic behaviour of £,r{t,r) at t — > -l-oo. If £,ruj{r*) is defined by an analytic 
continuation into the whole cut complex w-plane the contour of integration in (^) can be deformed down to 
lower values of Im lu by going around all the poles of ^rui {r* ) and alongside the branch cuts (see Fig. 2) . The 
part of the integral along the closing line Imw —00 vanishes for i > and the integral in ( |36| ) will be equal 
to the sum of residues of the poles at uj = LOnm and at cj = — fc and the contributions from both sides of each 
cut attached to the singular points cj = a>t, described above. Thus we obtain 

Ut, r)^iY, resU=.._ (^™e-*"*) + i res^^^fc (^^e"^'^*) + (46) 



-y 



+00 +00 
■I 
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where we have made the substitution lo — lui, = —ix and denoted the value of on the left side of the branch 
cut as £,1^, while on the right side as It is necessary to emphasize that the procedure of finding ujnm having 
Imcunm < depends on the arbitrary choice of the branch cuts, i.e. if someone do another cutting (which 
corresponds to another bypassing procedure when integrating in r-plane, say, to go around r^(a;) in the case A 
on Fig. 3), he will probably find out new eigenfrequences and lose the old ones. This indicates that we are not 
able to decompose the integral in (|3^) into the sum (^6|) by the only way: with the cutting being changed, the 
contribution from poles may converted partially into the contribution from cuts and vice versa. Formula ( ^ ) 
is valid for all t > 0. Because of the factor e~'"* rapidly decaying with t , asymptotically for t — > oo the main 
contribution to the £r{t) will come from the poles having maximal imaginary part and from real valued tOb- 
Terms coming to expression ( ^6| ) from the residues of the poles will be proportional to exp(— ia;„„iO 

Ut, Oldiscr = fc)e"^""""=)*, (47) 



and corresponds to the discrete spectrum of w. From formula ( |39| ) and expansion (45) of gi(l,a;) near the pole 
Unm it follows that the coefficients Enmif) in ( |47| ) must be proportional to gi{r). Hence, Enmif^k) are the 
eigenmodes for eigenfrequences LOnmik). They are the solutions of the edge problem for uniform equation (^) 
with the integration along the contour bypassing Alfven resonant points r — r^(a;) in complex r-plane as 
described above. Second term in the expression ( p6| ) is proportional to exp{ikt), it's radial profile is not an 
eigenmode and depends on the initial conditions. As it has been already mentioned, this term is the result of 
the degeneracy of one of the Alfven resonant points arisen due to the Bq^ — const. 

For real r on the interval (0, 1), which only are meaningful for physics, radial eigenmodes are continuous 
functions when Imwnm > 0, have logarithmic type singularity at r — r^(a;„m) when Ivaujnm — 0, and have a 
discontinuity at the point r = ryi(Reaj„m) of intersection of the branch cut with real r axis when Imaj„„i < 0. 
In the last case the position of discontinuity depends on how we do the cutting of r-plane (which is related to the 
cutting of oj-plane for each < r < 1), and is r = r^(Rea;„m) only by our particular agreement on the cutting 
of w-plane. Radial eigenfunctions are well determined only for eigenfrequences with nonnegative imaginary 
part, i.e. for unstable or neutrally stable modes. This is not surprising because, as has been already mentioned 
above, one can not even uniquely determine the whole multitude of stable eigenfrequences uJnm- Such situation 
always takes place in the problem of the stability of force-free relativistic jet, because, as seen from (^2|), Alfven 
resonant point lying in the interval (0, 1) can be always found at some real frequency uj for any k, m ^ 0, 
and for any choice of the function il^(r). All above treatment of the initial value problem involving Laplace 
transformation shows that this conclusion on eigenmodes have to be true for the stability problem of any ideal 
hydrodynamic flow with the equilibrium conditions depending only on one space variable (for the problem being 
reducible to the ordinary second order differential equation), whenever resonant surfaces of the perturbation 
with the flow or with the characteristic waves of the medium (sonic, Alfven, slow magnetosonic) do exist for real 
oj. Particularly, as pointed out in BIB, for cylindrical nonrclativistic MHD flow singularities with the Alfven and 
slow magnetosonic waves (when A = or S* = in the notations of BIB) exist only when ui is real. Moreover, 
these resonances are the logarithmic branch points for the solutions of radial differential equation. Therefore, all 
our consideration is directly applicable to that case, with the only correction for more complicated structure of 
the factor ahead of dS^r/df and dp^^/dr in BIB equations (3) analogous to our equations (p8|), and, hence, more 
complicated structure of analytical continuation of in the lower w half-plane due to the enhanced number of 
possibilities for singular points to merge together. In all previous investigations known to us the authors either 
restricted themselves to finding only unstable modes in their numerical computations (Torricelli-Ciamponi & 
Petrini 1990, Appl & Camenzind 1992) or chose some particular uniform profiles for vo(r), Bo(r), equilibrium 
pressure po(^) and density pQ(r) such that the factor ahead of dS,r/dr and dpi,/dr became independent on r, 
and all resonances disappeared (Turland & Scheuer 1976; Blandford & Pringle 1976; Hardee 1979; Cohn 1983; 
Payne & Cohn 1985). Therefore, it could be possible to integrate radial equation along the real r axis only, and 
in simple cases even obtain the dispersion relation I?(w, k, m) = expressed through well known mathematical 
functions. When finding stable modes it is necessary to bypass some of the singular points in the complex 
r-plane according to the procedure described above even if they are not on the real axis interval (0, 1). 

Consider now the third term in ( |4^ ) rising from integration along the branch cuts in w-plane. From (|39|), 
( ^0|) and (^) it follows that the main terms of the ^^(w) expansion in lo — lda in the vicinity of the branch 
point LOAif) are 



£.ru,ir) = C3o(wA)log 



"^{u - ujA{r)) 



, . , 50(1, wa) 1 , 

Clo(t^A) + —7^ 7 log 



r 



-(w - iOA{r)) 



(48) 



with the dropped terms giving nonzero contribution to (Uq) not greater than {uj — loa) Xo'giuj — uj a) ■ Substitution 
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of (|4^) into (^6|) after performing calculations gives us two leading terms of for t — > oo 

Ut, r) = [a,(r) log [y^^t) + «2(r)) + . . . . (49) 

At different points r the perturbation has different frequencies LOA{r) and different phase velocities 

This corresponds to continuous spectrum of w, since for fixed values k and m all frequencies w, with r^(w) 
lying between and 1, are exited. As in the cases of the sheared noncompressible hydrodynamic flow under 
the presence of the critical surface (Timofeev 1970, for a review) and diffuse MHD pinch (Kadomtzev, 1988) 
the perturbation as a whole can be treated as slitted onto a number of localized, singular in r perturbations 
propagating with the local Alfven velocity ( [50| ) . The decaying of perturbation observed from (|4|) is due to the 
phase mixing of neighbouring localized perturbation when they propagate with different velocities ([50|). Because 
for the flow considered in this paper continuous spectrum always present and is real the asymptotic for t oo 
behaviour of perturbations is determined either by an eigenfrequency with Re LVnm > (exponentially growing, 
unstable case) or by real eigenfrequencies if the former is absent (oscillations, neutrally stable case), or, at last, 
by continuous spectrum (E3) (decaying with time, stable case) if there is no nonnegative eigenfrequencies. 



3 Results of numerical computations 

To find eigenfrequences Wnm it is necessary to solve uniform equation (|3^) . Even in the simplest case fl^ = const 
this equation already has too many singularities to be solved by the well known hypergeometric functions and 
we are led to either doing numerical computations or finding conditions and small parameters under which it 
can be simplified. The results of numerical computations are presented in this section. 

In accordance with the properties of described above we adopted the following procedure in numerical 



computations. We integrated in r pair of equation (g8|) instead of ( |32| ) because they have simpler coefficients 
and do not contain additional singularity at C2 = 0. Integration was based upon five order Dormand-Prince 
method. For solving eigenvalue problem we applied shooting technique. For every k, m and chosen initial value 
of uj we started from the point very close to r = with the initial condition oc r'™'"^ and sought after 
the ^r-|r=i being equal to 0. If Imw > integration should be performed entirely along the real r axis. For 
Imo; negative, equal to zero or very small positive we passed round the zeroes of A in the complex r-plane in 
accordance with the consideration in section 2: it is necessary to find the position of zeroes of A in the r-plane 
for frequency Woo having real part the same as w, but the imagine part very large positive, then plot their paths 
in r-plane when diminishing Imw from -|-oo to the value at hand. If the path of a zero intersects with real 
r axis at a point in the interval (0, 1) than this zero is needed to be passed round, and we have to integrate 
equations ( p8| ) along the contour going around the rA(w) (see Figs. 4 for different cases of paths). In actual 
computations it occurred that the eigenfrequencies conm have small negative imaginary part (module of it is less 
than 0.1). If we consider only uj with |Ima;| ^ 1 than the points rA{uj), which are needed to be passed round, 
will lie close to the real r axis. The intersection point rA{uJo) of a path of rA{uj) with the real r axis, when 
changing Imw, will be located close to the rA{uj). Therefore, numerical procedure can be simplified because 
there is no necessity to follow full path of each r^(w) when Imuj runs from +00 to the value at hand. We 
actually did the following: took real ujq = Rew, found all rA{i^o) belonging to the interval (0, 1), then found all 
rA{(^) and pass round only those of them, which are the nearest to rA{uJo)- In doubtful cases we did the full 
procedure of finding which points from all multitude of r^(aj) must be passed round, being provided for the 
separate code. The value of at r = 1, which was being achieved to be 0, does not depend on the particular 
contour of integration, provided the points rA(w) are passed round properly. Hence we can suite the contour 
of integration for programming convenience. We chose the rectangular contour when going around the points 
rA^oj), and, in order to save the accuracy of computations in the vicinity of singularity r — ta, used to passed 
round also the points rA(w) when Imuj was very small positive value. 

We see that for small Imujnm the rule for bypassing singular points in equations ( |2^ ) coincide with that 
used in calculations of Landau damping of waves in plasma medium when doing integration in velocity space 
(Lifshitz & Pitaevskii, 1979). There are two differences, however. First, in our case integration is performed in 
the interval (0, 1), while in the case of Landau damping along the whole real r axis from —00 to +00. Second, 
not all singular points of ^rui in ^ plane are poles, but the branch points as well. As described in section 2 these 
lead to the existence of continuous spectrum for the problem considered in present work. 

It occurred that lmujnm{k) < for all fc, ni, and 3 functions of il^{r) involved in calculations, so the jet 
is stable with respect to helical perturbations as well as with respect to axisymmetric perturbations. Actually, 
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computations were performed for — 10 < fc < 10, —3 < to < 3, and for the first 3 radial modes for each k and 
TO. The dependencies 

n^ir) = n{l~r^), n^{r) ^ne-''cos(—r) , and 



(r) = n 



2 



a + b , 1^ ^ 



ah 

2 3 



where 17, a, b are constants, < a < 0.5, 0.5 < 6 < 1, were tried out. Constant fl in the expressions for fl^ was 
ranged in the interval from 0.1 to 20. 

Our model does not itself provide us any information about the function fl^{r) except that = 0. The 

angular rotational velocity of magnetic field lines has to be found from the consideration of the jet origin. We 
adhere the viewpoint that the part of the jet closest to the symmetry axis is connected by magnetic field lines 
directly to the black hole and can be described in the frame of force-free approximation (jets emerging from 
the accretion disks seem to be not the force-free, because of the kinetic energy of the mass flow is comparable 
to the Pointing flux, Pelletier & Pudritz (1992)). Therefore, to determine U^(r) one should solve the equations 
governing stationary two dimensional axisymmetric structure of magnetic fields in the vicinity of a rotating 
black hole taking into account the unavoidable process of particle creation there (Blandford & Znajek 1977, 
Takahashi et al. 1990, Nitta, Takahashi & Tomimatsu 1991, Beskin, Istomin & Pariev 1992(b)). This is stiU 
unresolved problem. Simple model describing the force-free magnetic field in the vicinity of a slowly rotating 
black hole was developed disregarding the effects of e"'e~-pair creation in Beskin, Istomin & Pariev, 1992(a). 
The dependence 

f7^(r)=fi-^^i^= (51) 
1 + Vl — 7" 

with unspecified constant H, follows from that model. To be closer to reality we adopt for computations of 
disturbances propagation along the jet the function 

n^{r) = n{l-r^), (52) 



which resembles the function ( |5l| ) for < r < 1, but is the entire function in complex r-plane, contrary to ( |51 

Following the procedure outlined we have calculated the dispersion curves LUnm — ^nm{k). First three 
branches (n=0,l,2) of them for Vl^ given by (|5^) and to = 2 are shown on fig. 5 and fig. 6. On fig. 7 we show 
the dependence of the real part of u!nm{k) for m = —1. In this case, because of the absence Alfven resonance 
surface inside the jet, imaginary part of ujnm{k) is always equal to (continuum spectrum does present, of 
course). If some 3 values fc, w, and to are a solution of the eigenvalue problem than the values — fc,— w*, and 
— m will be a solution too but for the complex conjugated function so we depicted only the branches of 
ujnm = w„m(fc) having Rew > 0. Those having Rew < can be obtained by the reflection of fig. 5 and fig. 7 
with respect to the coordinates origin. On fig. 8(a) we plotted an example of radial eigenmode when there 
is no Alfven resonant point on the interval < r < 1, on fig. 8(b) the same for the case when resonant point 
exists. Note the localized character of the mode in the last case. In order to find physical reason for decaying of 
the discrete modes having Alfven resonant point we calculated Pointing flux S = l/47r(E x B) for each mode. 
On fig. 8(c) the radial component of energy fiux rS'r, calculated for mode shown on fig. 8(a), is depicted as a 
dependence from r. It is seen, that Sr < 0, i.e. the energy flows to the jet axis, and, which is more important, 
Alfven resonant surface is a drain of electromagnetic energy of the mode. This energy is converted into the 
energy of the mean magnetic and electric fields, so one should expect that near Alfven resonant surface strong 
amplification of the magnetic and electric fields does occur. However, the consideration of this process is based 
on the second, nonlinear approximation and is beyond the scope of the present work. 



4 Boundary layer analysis 

In this section we consider the spectrum of eigenfrequencies ujnm near the merging point Uc of two VAii-^)- At 
this point A(r,io) has at least quadratic zero in r. We assume that A (wc) ^ and the point rA{i-Oc) does not 
coincide with and 1 (prime denotes partial differentiation with respect to r at the point r — Vc). To find 
points LOc{k,m) it is necessary to solve the couple of equations A(rc,i-Jc) — 0, A'{rc,i-Uc) = 0. The solution can 
be represented in implicit form as 

fc = -mn^ir^) + -^^^lil^il + n^'{n)rl), (53) 
= mn^in) + ^^^li^il - n^\r^)rl). (54) 
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There are two cases of the sohition of the system (|5^), ( ]54| ) for each k and m. In the first case, when the 
solution Tc of ( p3|) is real, ojc determined by ( ^4[ ) will be real. In the second case, when (js^) has a pair of 
complex conjugated solutions Tc and r*, equation ( |5^ ) will give us a pair of complex conjugated values Uc and 
Lu*. As described in Subsection 2.3, contour of integration of equation (32) can never being clamped between 
merging points r/i(w) when Imwc > 0, so we are able to move it away out of the neighbourhood of u>c- Hence, 
uJc with Imwc > is regular point of and can not be an accumulation point of the poles. Perturbations 
with luiLU < are exponentially decaying with time and do not contribute to the asymptotic of S,r{t) at i — > oo. 
Consequently, the spectrum Unm near the real points coc is of primary importance for the stability problem. If 



= 0. 



u)c and Tc are real, the frequency uic coincides with the edge of the Alfven continuum, i.e. '^^^^^ 

Below we perform boundary layer analysis of equation (|3^) in the vicinity of Vc similar to that was done 
in BIB for the edge of the slow wave continuum in the nonrelativistic MHD stability problem of nonrotating 
cylindrical flow. But in contrast to BIB, who always integrated radial eigenvalue equation along real r axis, we 
adhere the rules of contour deformation in the complex r-plane, which are described in Section 2. At the point 
r = Vc, the eigenvalue problem ( p^ ) has a singularity such that to lowest order in x ^ r — rc 

A~A"x^/2, 

Ci, C2, C5 ~ const, C4 ~ ~A"{uj + k)x^/2, 
(Ci/rCa)' ~ const. 



Thus, to lowest order in x, equation ( |3^ ) becomes 

d 



d_ 

dx 



where 



dx 



A" 



C5 + rC2 



rC2 



Setting r^r oc x'^ , we obtain for the characteristic exponent 



-1/2 ± ^1/4- D. 



(55) 



(56) 



(57) 



If 



D > 1/4, (58) 

than will have infinite number of oscillations in the vicinity of r = Tc and j^^l oc r~^/^. If Z) < 1/4, 
than two linear independent solutions will be of power type without oscillations. As we will see below, in 
the case 13 > 1/4 there exist an infinite sequence of eigenfrequencies ujnm- Asymptotically, for 71 ^ 00, the 
eigenfrequencies converge geometrically toward the marginal point uj = ujc- The eigenfunctions corresponding 
to these eigenfrequencies have the higher number of oscillations in the vicinity of r = Vc, the higher is the 
number n, and the limiting solution for has an infinite number of oscillations. In the other case, D < 1/4, the 
picture outlined fails to be true, and the lowest-order boundary analysis can not provide any information on 
eigenfrequencies. Using definition (|5|) and relations (p3|), ( [54| ) oscillations condition ( p8| ) can be transformed 
to the form 

d_ 

dr 





d 


" 8 




dr 


r 

(f^^^r2)' 



< 0, 



(59) 



where all derivatives are evaluated at the point r = Vc- We see that this condition depends only on the position 
of the point rc, and does not depend on k and m separately. 

Let us now consider the frequency lu slightly different from ljc 

uj ^LUc + A, (60) 

where |A| ^ 1. In the vicinity of a; = Wc, r = Vc the main terms of expansion A in uj and in r are 



A 



(1 + f7^"r2)A + 2mn^" - {iVc + k){n''^r^)"] x^/2 



(61) 



Note, that the equation ( |55| ) is invariant under rescaling of x. Furthermore, from (|61|) it follows, that the 
rescaling of x leaves the lowest-order in A and in x approximation of radial mode equation ( p2| ) invariant if 
A is scaled the same as x^. Therefore, boundary layer of thickness x oc A^/^ arises near the point r = rc in 
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the equation (p^). To obtain "inner layer" equation we introduce inner rescaled variable X by the following 
definition 

1/2 




1 2 



-^ = — (-7T- . ^.2 . I , (62) 



where we determine v A such that the —tt/2 < argv A < 7r/2 and accept for the argument of the square root 
from the expression enclosed by parentheses the value 0, if this expression is positive, and the value tt/2, if this 
expression is negative. With such definition the lower-order inner-layer equation becomes 



d 

dX 



D{r^r) = 0. (63) 



This is Legendre's equation with singular points X — 1 and X — —1. At these points r^r has logarithmic 
singularities, as they are merely two closely spaced first -order zeroes of A in r. Equation (|6^) shows that within 
the "inner layer", \X\ ^ 1, the finite frequency shift is essential. Furthermore, for \X\ ^ 1 and |a:| ^ 1, 
there exists an "intermediate" region where the frequency shift becomes negligible, and where the inner-layer 
equation ( p3[ ) approaches the small |a;| limit of the external equation (^5|). The boundary layer analysis 
consists of matching the solutions of the external equation ( |55| ) for \x\ small with those of the complete inner- 
layer equation ( |63|) for \X\ large, using (|6|) to relate |a;| and \X\, to obtain an equation for A. Figs. 9 (a), (b) 
show the complex x and X planes for the case when ImA < 0. On Fig. 9 (a) mil^ /D < 0, on Fig. 9 (b) 
mft^ /D > 0. For definiteness, we assume that the contour of integration of equation ( ^2|) approaches the two 
merging points ^^1(0;) along the real x axis from x < and moves away them also along the real x axis to a: > 0. 
Let us introduce the notation for characteristic exponent (|57|) 



i/ = -1/2 ibis, ^JD- 1/4. (64) 

Parameter s is real for oscillatory solutions and imagine if the otherwise. Then, the solution of the external 
equation ( ^5| ) in the vicinity of r = is for a; < 

r^r = a+x-i/2+- _,. a_a;-i/2-^^ (65) 

where a-^^/a^ is fixed by the external solution in order that the boundary condition at r = be satisfied. 
Similarly, for a; > 0, 

r^r = b+x-^^^+'' + b-x-^/^-'', (66) 

and b^/b^ is fixed by the external solution in order that the boundary condition at r = 1 be satisfied. Matching 
the solution of the internal equation with the solution (p5|), ( |6^ ) of the external equation leads to asymptotical, 
when \X\ —^ 00, expressions for inner-layer solution 



r-T- N-l/2+'iS / r-^ s -1/2-lS 

V AX j + 6_ (VAX) when ReX +00, 

V AX ) + a_ f V AX ) when ReX ^ -00. 



Here we denote 



~ ,12 mn'" 
A = -A- 



Since the scaled inner-layer equation ( p3[ ) is independent of A, so is the connection matrix between the expansion 
coefficients for r^r in the two asymptotic regions ReX +00 and ReX —^ —00 divided from each other by 
branch cuts attached to the points X = +1 and X = —1 (if the contour of integration goes between these 
points). Thus, 




V21 V22 



(67) 



Here Vij only depend on D. The explicit form of vn, V12, V21, V22 follows from the solution of the inner-layer 
equation (p3|). The general solution of this equation is 



r^r = AP,{X)+BQ,{X), 
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where P,^{X) and Q^{X) are Legendre's functions of the first and second order respectively with v from (p^). 
The main terms of the asymptotic expansion of Pi,{X) and Qi,{X), when \X\ — > oo, are (Abramowitz & Stegun, 
1970) 



UX) 



0Fr(i/2-is) 

-T{l/2 + is) 

' r(i + zs) 



^(l/2 + ^s) 



(2X) 



-l/2-is 



The form of these expansions just match the expansion (|6^), (|66|) of the solution of the outer-layer equation. 
Therefore, considering branch cuts as shown on Figs. 9, involving the relations between Py and Qi, on the 
different sides of a branch cut, and properly choosing the branches of functions x^^l'^'^^^ in the asymptotic 
expressions, one can obtain the coefficients of connection matrix (if the contour of integration goes between the 
points X = +1 and X = — \) 



Wll 



= -e ^'"*(1 + cothTTs), W22 = e^'"*(coth7rs - 1), Wi2 = %/7rse^ 



r(l/2 + is) 



V2\ 



se 



-2irs 



TT \r{i/2~ is) 



Then, taking into account that ~ b^/b_ and i?_ = a^/a_ are fixed by boundary conditions, we find that in 
order for the system ( |67[ ) to have nontrivial solution for a+ , a_ , 6+ , 6_ , the following dispersion relation must 
be fulfilled 



This is a quadratic equation on A'* with the solutions 



± 



y/R+R- sinhTTs ^ \^ sinh 



'R- 



^cosh TTs 
sinh^ TTS 



Because of the multivaluedness of the raising to a complex power there exist two infinite sequences of the 
solutions A 

1 . 27rn i 
-argC± log|C±l 



A, 



"(±) 



exp 



(68) 



where n is an integer number. If s is real (s > 0) than with increasing n the module of A„(-|-) will be decreasing 
and the sequences of eigenfrequencies uJnm = ij^c + A„(-|-) will converge to marginal point ujc- Note, that the 
analysis in this section is only valid for small A, i.e. the expression (Em is accurate only when n is large and 



positive. If s is imagine and oscillations condition (58) breaks down, this will not be the case. If D < 1/4 
than An{+) a-nd A„(_) from (^8|) have the same modules for every n, with only the argument changing with n. 
Hence, these eigenfrequencies are improper, and in reality there are no eigenfrequencies in the vicinity of Wc- 

Without knowledge about i?+ and i?_ one could say nothing about the sequences of A„(-|-) except their 
existence. In turn, the values of R+ and i?_ depend on the solution of full equation (^2|) with possibly existing 
another points r^(w) needed to be gone round when integrating (^. 

At last, we mention that in the case when contour of integration does not go between the points X — +1 
and X ~ —1, but go round the whole couple, than the coefficients of the connection matrix in ( |67|) become 
Vii = V22 — Ij vi2 = V21 — 0. Hence, one can no longer obtain any information on A from (|67|), but natural 



conclusion that 64 



and 6_ 



follows due to that both branch cuts are on the same side from the 



contour of integration. Consequently, as already pointed out in the beginning of this section, such points ojc are 
not remarkable. 



5 Possible astrophysical implications 

We believe that strongly magnetized inner parts of the jets can be directly connected to the central supermassive 
black hole by magnetic field lines. Let us estimate wether such jet can be force-free. We assume the black 
hole have typical values of its mass and rotating parameter: M ~ 10^ Mq, a = J/M < 1. On the distances 
from black hole comparable to the Schwarzschild radius strong magnetic field is expected to be of the order 
of 1Q*G (Begelman, Blandford & Rees, 1984). Than the electron density of ric ~ Icm'^^ is enough to screen 
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the longitudinal (along the magnetic field) electric field component so that the MHD approximation becomes 
possible. In the inner part of the flow connected with the black hole by magnetic field lines the particle density 
n cannot much exceed the value ric because the particles constrained to move along the magnetic surfaces 
do not escape the black hole and the only source for replenishing them is e+e" pairs production. However, 
particles creation in the black hole magnetosphere seems to be possible only in the presence of the longitudinal 
electric field, which vanishes for n » ric- Therefore, an equilibrium between the particles outflow fed the jet 
and their creation is established for particle densities of the order of ric- The necessity of particles creation in 
the magnetosphere of supermassive black hole was first pointed out by Blandford & Znajek (1977). In this case 
the energy density of the magnetic and electric fields is 10^^ times greater than the rest energy density of the 
e^e~ pairs. That makes the force-free approximation adequate in the vicinity of supermassive black hole. 

Apparently, the force-free approximation is also valid for the inner parts of the jet, which are close to the 
axis of symmetry and are connected with the black hole. To show this consider the conservation of the particles 
flow and the conservation of the current in the process of possible recoUimation of this inner part of the jet 
from the sizes of the order of the black hole event horizon to much greater. Continuity of the particles flow 
implies that 7n oc l/B?, where R is the radius of the jet, 7 is the Lorentz factor of the plasma flow, n is the 
particles concentration in the reference frame comoving with the plasma flow and the velocity of the particles 
is relativistic everywhere. Conservation of the current flowing inside the jet leads to that the magnetic fleld 
after recoUimation will be predominantly toroidal and scale as B oc \/R. We see therefore that the ratio 
B'^ / ATrmnj'^ oc I/7, so after recoUimation to larger radii (say, parsccs) the jet remains to be force-free unless 
the Lorentz factor of the flow will not reach an unbelievably high value lO^''. Therefore, we hope to apply our 
consideration for those jets, whic;h are believed to be electron-positron, and are likely to be force-free. 

The remarkable feature of the dispersion curves Reuj{k) is that they have a minimum at some k = kmin and 
kmin 7^ (see Figs. 5 and 7). At the same time waves damping lmuj{k,miri,) is either small (it never exceeded 
0.1 in our computations) or even equal to for modes with m < 0. Because of these, the perturbation with 
k « kjnin do not propagate since the group velocity duj/dk vanishes for k = kmin- In contrast to the waves 
having k ^ kmin this wave packet undergoes only diffuse broadening due to the finite value of cPuj/dk"^ for 
k = kmin- It means that such oscillations form the "standing wave" with the wave vector kmin- Such "standing 
wave" does not transmit any energy because of vanishing group velocity. We use quotation marks to name this 
phenomenon in order not to confuse it with the well known in many fields of physics standing waves, which are 
the linear superposition of two progressive waves. The amplitude of the "standing wave" will be larger than 
the amplitudes of other waves because it experiences a dispersion spreading only. This phenomena is caused by 
the fact that the oscillations with wave vectors less and greater then kmin propagate in the opposite directions. 
The phenomenon of "standing wave" takes place for axisymmetric perturbations as well (IP). 

After a long time after initial excitation the pattern of disturbance is formed with the wave crests moving 
with the phase velocity uj{kmin)/kmin- In our numerical calculations this velocity was always greater than the 
velocity of light. A "standing wave" perturbation occupies progressively growing with time as ex t^^"^ part of 
the jet. The edges of this pattern move with velocity oc which is less than the velocity of light. At the 

same time amplitude of the perturbation inside the "standing wave" region is decaying with time as oc 
It is impossible to transmit any information by moving wave crests faster than the speed of light. If one has 
relativistic electrons emitting synchrotron radiation inside the magnetic configuration dealing with (which is the 
case for extragalactic jet), than this pattern will be visible. This provide us with the new type of superluminal 
source. Now according to well known formula one can calculate the observable velocity of such superluminal 
source in the projection onto the plane of the sky 

vsinO 

Kbs = 1 ^, 69 

1 — v/ccosO 

where is the angle between the jet axis and the line of sight of the observer, v is the velocity of the superluminal 
source. In Fig. 10 the dependence of V^bs from 6 is depicted. If ^ < ^* = arccosc/u than Vohs < 0, i.e. the 
apparent motion of knots will be reversal, in the direction opposite to the wave vector kjj^jjj. Observer will see 
superluminal motion (| Vohs |> c) if 

c " 
arccos— ;= 45° < 135° — arcsin 



\/2u ^/2v' 

Note, that the velocity, which enters into the effect of relativistic beaming, is the velocity of matter flow inside 
the jet rather than the velocity v. Consequently, the brightness of superluminally moving knots will not be 
affected by their motion. 

The source of perturbations may be instability of accretion disk around the central black hole, which affects 
the magnetic field in the vicinity of a black hole, and, therefore, lead to the excitation of the disturbances at 
the base of the jet. If the jet is oriented close to the line of sight {0 < 0*), the observer will see chain of knots 
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moving backward to the core. In the counterjet, if observable despite the fading due to the effect of relativistic 
beaming, knots will move outward from the core, as seen from Fig. 10. Another possibility is the excitation of 
perturbations in the jet at the hot spot, where the jet is ended. Though we consider only stationary equilibrium 
jet structure, the velocity of the advancing of the jet into the extragalactic medium is believed to be only mildly 
relativistic 0.2 c) (see, e.g. Begelman, Blandford & Rees, 1984), so one can admit that the phenomenon of 
"standing wave" may take place in that region as well. Here knots in the jet, provided that 6 < 9* , will move 
from the hot spot in the direction of the core, while in the counterjet they will move to the hot spot out of the 
core. If 6* > 6** then the observable direction of motion for both jet and counterjet coincides with the projection 
of the real motion of knots onto the plane of sky, i.e. outward from the core and outward from the hot spot. 

From expression (|6^) it follows that ii ^ 0* , then |T4bs| oo having different signs for < 6* and > 0* . 
The equality of V^bs to infinity means that the whole jet splashes simultaneously throughout all it's length. In 
reality this can not come true. The jet radius changes along it's length, jet are usually slightly bent. These 
two reasons limit the value jVobsl by some high but finite value. If the viewing angle of the jet is near 0* and 
changes along the jet such that there exist pieces of the jet viewed at an angle less than 0* and greater than 
0* , then pairs of knots can be observed which either collide and vanish or emerge and move in the opposite 
directions, depending on whether changes along the direction of real motion of wave crests from being greater 
than 0* to being less than 0* (in the case of merging knots) or from being less than 0* to greater than 0* (in the 
case of newly born pairs of knots). When the wave crest of a "standing wave" pattern passes that piece of the 
slightly bent jet, where Ri 0* , observer will see that the knot, corresponding to that wave crest, is stretched 
along the jet and it's total luminosity (not the brightness) becomes considerably larger. Then the next wave 
crest passes through that place of the jet, and observer sees the next flash. This can be the reason for the 
quasi periodical splashing of the innermost parts of the jets with modest viewing angles « 0* of the order of 
45°, while usually such bursts are explained by relativistic beaming of the moving knots along slightly curved 
trajectory in the jet having no more than a few degrees (see recent paper by Camenzind & Krockenberger, 
1992). However, careful examination of these effects needs considering physical equation describing disturbances 
propagation inside curved jet for finding dependence v{0), not only kinematic picture. This is beyond the scope 
of the present paper. 

It is the task for radioastronomers to detect such "natural" (not due to the effect of projection) superluminal 
motions. Baath (1992) reported about three epoch observation of one component in 3C345 moving inward to 
the core, but he writes that the significance of this observation still remains to be verified. Hardee (1990) 
proposed another scenario which can lead to observation of backward motions of the intersection points of the 
shocks in nonmagnetized jet. In the frame of our model periodical structures moving backward to the core may 
be observable while Hardee's model predicts an isolated knots. 

6 Summary 

A considerable amount of extragalactic jets are extremely well collimated and extends over the distances tens 
times longer than their radii. Bright well known example of such jet is one emerging from the galaxy M87. 
The problem of extraordinary stability is long standing problem. In this work we have shown numerically that 
a jet with a longitudinal current is stable within the force-free approximation for all velocities of longitudinal 
motion and for wide range of the velocities of rotation. We consider initial value problem for linearized set 
of relativistic equations describing disturbances propagation inside cylindrical jet. By using Laplace transfor- 
mation we find asymptotical behaviour of perturbations over long time since initial excitation. The stability 
problem is reduced to eigenvalue problem for radial modes. It was shown that there exist Alfven continuous 
spectrum of eigenfrequencies and discrete spectrum having the accumulation point where two Alfven resonances 
coincide. Numerical calculation shows that all eigenfrequencies have been computed have negative or equal to 
zero imaginary part. This means the stability of the jet with respect to helicoidal perturbations as well as to 
axially symmetrical or pinch ones. The physical reason for stability is that there is a shear of the magnetic 
field because of changing the curling of the magnetic field lines with the radius (the absence of the shear would 
imply ri^(r) cx l/r, which leads either to the value of rotational velocity il^r + KBq^ being undefined on the 
symmetry axis of the jet, if 1 + KBq^ ^ 0, or flow velocity of the matter being equal to the speed of light on 
the symmetry axis, if 1 + KBq^ = 0, both cases have no physical sense). Because of the fluid pressure is low 
compared to that of the electromagnetic held, even a small shear stabilizes the motion perpendicular to the 
magnetic held lines and prevent the development of instability. 

We also find that the dispersion curves uj = Lj„„i(fc) have minima for certain values of fc = k^in- This means 
that such oscillations form a "standing wave" with a wave vector fcmin. The wave crests of this "standing wave" 
are spirals moving along the jet with the velocity exceeding the speed of light. The amplitude of the "standing 
wave" will be larger than the amplitudes of other waves because it experiences a dispersion spreading only. 
An example, illustrating this behaviour of perturbation, was presented in IP. Relativistic particles emitting 
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synchrotron radiation in the jet make the perturbed magnetic field visible for the observer. Thus we obtain a 
new type of superluminal sources, having their physical rather than only observable velocity superluminal. The 
observation of the chains of knots moving backward to the core of the extragalactic radio source will be the 
evidence that the phenomenon described does take place in Nature. The phenomenon of "standing wave" may 
also account for the quasi periodical flashes of the central compact sources. 
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Appendix A 



Here we list the expression for different components of the perturbations of the magnetic field, electric field 
and velocity field in terms of Bri and ^r- Force- free approximation and Bzq = const are adopted. We use the 
abbreviations (p2[)-(24) of the main text. For magnetic field 
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For electric field 

Eri — 



Bzi 



i uj + k 
S F 



dBri 
dr 



{u! — k — mU,^) — Bri- 



dr 



F 



1 



Bri{u! — k + mfl^) 



i dBri 1 
'S dr F 

{-koj + k^ - 



-UJ h rk{u + k) + Vt^ 

r r 



toI^-^w + nC jr^) + B 



moj 



dn^ r{uj + k) 
[ -2u? -ujk + k^ 



Ezl — Brifl 



UJ 



k + mil^ 

UJ + k 



m 



(A2) 



(A3) 



(A4) 



(A5) 



In the approximation of ideal force-free plasma the component of the velocity parallel to the magnetic field line 
does not enter into the basic equations (|l|), (||) of the main text and can be free. Therefore the component of 
the velocity first order perturbation parallel to zero order magnetic field remains also to be free. Another two 
components are 
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To obtain expressions for all components of the electric and magnetic fields and velocity we substitute in 



(A1)-(A7) for Bri it's value from equation ( |19| ) of the main text, which express Bri via S,r- Thus, we lead to 
the following 

Bri = iBzoF^^r, 



B,^ = -Bzo\^^r^+&-^^{n^r) + ^ 



Bzi = B 



dr r ) rS 



Eri=Bzoin^r^+ir^in^r) 



E^^ = -iBzQ{uj - ■mQ^)^r, 

Ezi=iB,o^^r{u; + k)^r, 
Vri = —iBzoiio — mfi^ — KBzoF)^r, 
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(A9) 
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(All) 

(A12) 
(A13) 
(A14) 



dn 



Bzqv^^ - B^^vzi = -Bzo\ ''^r-j^{l + KBzo) + -^{mQ^ - lu + KBzqF) 



(A15) 



The function £,r{r) has logarithmic singularity at the point r = r^, which is the simple zero of A (the case of a 
second order zero is considered in section 4), and is regular everywhere in the remaining part of complex plane 
r, except infinite point. Therefore, one can conclude from equations ( [A8| )-( A15) the following behaviour for 
disturbances of different physical quantities near Alfven resonant point r — rA {x — {r — rA)/rA) 

Brioclogx, cx -, Bzi(x-; 

X X 

Erioc-, cx logs, £^ziCxlogx; 

X 

t;riOclogX, BzqV^^- B^^Vzi (X -. 



It might appear from equations (A£)-( All ), ( A15 ) that the quantities B^f,^, B^i, E^i, B^qV^^ — iJ^gfzi are 
singular at the r — rs, which is the zero of S. However, this is not the case. Being regular at r = function 
can be expanded into power series in the neighbourhood of r = rg. By expanding second order differential 
equation ( ^ ) on one can find that t he co m binat ion rdS^r /dr{n^ {oj + k) — m/r^) — £_r{^^ {ijJ + fc) + m/r^)^ 
which enters into all expressions ( A9)-( All ), ( A15 ), becomes equal to when S — 0. Therefore, all physical 
quantities regular there. The only singularity may arise when A = Q, i.e. for modes from Alfven continuum. 

Appendix B 

In the general case, when Bzq 7^ const, one can perform calculations analogous to that was done in section 2 

when deriving the radial eigenmode equations (psl). The coefficient by the derivatives --r{T£,r) and is of 

' — ' r dr dr 

special interest, because it contains the singular points of the system of two first order differential equations on 

the displacement and the disturbance of the total pressure. It occurred to be the following 

A = BzliuJ - mn^f + (luB^^ + n^rkBzof - (^B^q + fcS.o)' . (Bl) 



Zeroes of A{ujAif) — give us Alfven resonant frequencies uja at a given value r. The expression (Bl) is 
quadratic in lj, therefore, for each r the equation ^ = has two solutions for ujA{r). If we suppose Bzq = const 
and use for B^^ the expression ^ with the sign '+', the equation (Bl) will become 



A = BzliuJ + fc) L - fc - 2TOf7^ + n'^^r^iuj + fc)j = BzlioJ + k)A. 



(B2) 



So A factorizes into two coefficient, one of which do not contain the r-dependence. These two coefficient are 
present in system (|2^). 
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Next, let us show that at any real r the quadratic equation ^ = has only real solutions uoa. The 
discriminator of this equation is 

V = i-mn^ + n^rkrf - (1 + T'')n^\^ + + (1 + r2) (^r + fc)' , 

where r — B^^/Bzq- The solutions uja will be complex if and only if I? < 0. The condition P < transforms 
to 

which, in turn, means that E'ro > Bo^. In stationary ideal MHD configuration electric field can never exceed 
magnetic field. This would imply the velocity of the fluid v being greater that the velocity of light. Naturally, 
it can be easily shown that any solution of the equation (|^), governing the stationary jet configuration, satisfy 
the requirements \E\ < \B\ (see IP). Therefore, V must be nonnegative value and the solutions loa are always 
real. Thus, in the general case of force-free cylindrical equilibrium Alfven continua are always real and do not 
lead to instability. 
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Figure captions 

Figure 1. Schematic representation of the equilibrium stationary configuration of a jet with a uniform poloidal 
magnetic field Bz- The frequency of rotation in the dimensionless units described in the beginning of Subec- 
tion 2.1 is fl^ = 10(1 — r^). The jet boundary for r = 1 and three magnetic tubes for r = 1/4, 2/3 and 9/10 
are shown. The magnetic field lines are spiralling on a magnetic tube. Since ri^(l) = 0, the total current 
through the jet is equal to zero and the magnetic field is purely poloidal both at the boundary and at the axis 
of symmetry. The curling of magnetic field lines is maximum for r = 1/ \/3, decreasing for smaller and larger 
radii. The density of the poloidal current jz is negative when r < 1/V2 and positive when 1 > r > 1/^/2. The 
electric field E induced by jet rotation is radial. The plasma velocity v along the magnetic tube consists of 
two components: rotation with angular velocity n^{r), and motion along the magnetic field lines with a speed 
V|| = K'B. We see that the rotation velocity rH,^ can exceed the speed of light c ( in our case the maximum 
value of rO^ is 20/3-\/3c at r = l/\/3 ); nevertheless the quantity v is restricted by c due to the existence of a 
predominantly toroidal magnetic field. The dispersion curves uj = uJnm{k) for perturbations of this equilibrium 
state are plotted in the Figs. 5,6,7. 

Figure 2. Schematic picture, illustrating analytical properties of ^rtj hi the complex w-plane. Sec explanation 
in the text. Brunch cuts are shown by thing solid line, the contour of integration shifted into lower ui half plane 
breaks up into the circles around the poles u = LOnm{k) and the part going round the branch points and along 
the attached branch cuts. This contour is shown by bold solid line with arrows. 

Figure 3. Schematic picture of the complex r-plane. The trajectories of r^(a') when changing ui along the 
paths shown in Fig. 2, are plotted by dashed lines, deformed contours of integration in cases C and B are plotted 
by bold solid lines, branch cuts (one for each case A, B, C or D) are shown by thin solid lines. 
Figure 4. The thin curves with arrows show the trajectories of 6 points rA(w) when changing Imw from +oo 
to some large negative value. The endpoints are marked by crosses. In all figures = 10(1 — r^), m ~ 1, 
A; = — 3. Figures (a), (b), (c), (d) and (e) correspond to the continuation paths A, B, C, D, and E in Fig. 2 
respectively. uja{0) = 17, ^^(l) = —3, two complex conjugated u>c are iVc = 3.3117 ± 0.1564i, r* is posed at 
0.5. In figure (a) Rcu = —5, in (b) Rew = 5, in (c) Rew = 3.305 , in (e) Reo; = 3.315, in (d) Rew = 17.5. On 
figure (d) 4 lateral loops are enlarged, because their dimensions are less than 0.01 in reality. 
Figure 5. The dependence of the real part of cOnm{k)- m = 2, fl^ = 10(1 — r^). Three lowermost branches of 
the dispersion relations are shown. They are distinguished from each other by an appropriate number. Straight 
lines are oj = k and oj = —k. 

Figure 6. The dependence of the imaginary part of a;„„i(fc). m = 2, = 10(1 — r^). Three curves correspond 
to those shown on Fig. 5 and are indicated by appropriate numbers. 

Figure 7. The dependence of the real part of ijJnmik). m = —1, = 10(1 — r^). Three lowermost branches of 
the dispersion relations are shown. They are distinguished from each other by an appropriate number. Straight 
lines are oj — k and uj = —k. 

Figure 8. Radial dependence of the second mode for (a) m = —1, k = —0.4, = 10(1 — r^), wi-i = 6.86 
and (b) rn = 2, A; = 1, = 10(1 — r^), wi2 = 8.54 — 5 • W~^i. The solid curve shows Re^^ and the dashed 
line Im^r- The solution was started so that is real as r — > 0. Unit of is arbitrary. Fig. (c) shows radial 

dependence of the radial component of energy flux, calculated for mode (a). 

Figure 9. Complex x and X plane, shown for the case argA = —tt/'3. In Fig. (a) /D < 0, in Fig. (b) 

mfl^' /D > 0. Contour of integration is plotted by bold solid line, branch cuts are shown by thin solid lines. 
Figure 10. The dependence of the observable velocity of the standing wave pattern V^bs on the angle between 
the jet axis and the line of sight of the observer. Phase velocity of the perturbation v is chosen equal to 3c. If 
the jet is pointed directly to the observer than ^ = 0, if it is pointed directly from the observer than 6 = 180°. 
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